function [ U, S, V, out ] = ActiveSubspace( D, lambda, para)

tol = para.tol;
maxR = para.maxR;
if(~isfield(para, 'initR'))
    initR = ceil(0.01*min(size(D)));
else
    initR = para.initR;
end

[m, n] = size(D);
[row, col, data] = find(D);

R = randn(n, initR);
R = powerMethod( D, R, 3, 1e-6);
[U, S, V] = svd(R'*D, 'econ');
U = R*(U*S);
B = V;
subtol = norm(data, 2);

spa = sparse(row, col, data);
spa(m, n) = 0;

maxIter = 2000;
obj = zeros(maxIter, 1);
RMSE = zeros(maxIter, 1);
Time = zeros(maxIter, 1);
t = tic;
for i = 1:maxIter   
    Z = partXY(U', V', row, col, length(col));
    Z = data - Z';
    spa = setSval(spa, Z, length(col));
    
    % approximate SVT
    Q = powerMethodMatComp(U, V, spa, B, 3, 1e-6);
    Z = (Q'*U)*V' + Q'*spa;
    [Ua, Z, Va] = svd(Z, 'econ');
    Sa = diag(Z) - lambda;
    Sa = sum(Sa > 0);
    Ua = Q*Ua(:,1:Sa);
    Va = Va(:, 1:Sa);
    
    % makeup active subspace
    B = filterBase(Va, V, 1e-3/i);
    [B, ~] = qr(B, 0);
    B = B(:, 1:min(size(B,2), maxR));
    A = [Ua, U];
    [A, ~] = qr(A, 0);
    A = A(:, 1:size(B, 2));
    
    % solve sub-problem
    Zi = (A'*U)*(B'*V)';
    [U, S, V, ~, subIter] = actSSsubproblem(A, B, lambda, data, row, col, ...
        subtol, Zi);
    obj(i) = subIter(end);
    subIter = length(subIter);

    % update X
    U = A*(U*S);
    V = B*V;
    
    % check convergence
    if(i == 1)
        delta = inf;
    else
        delta = obj(i - 1) - obj(i);
    end
    subtol = abs(delta)/2.5;
    
    fprintf('iter:%d obj:%.4d (%.2d) rank:(in %d out %d) inner iter:%d\n', ...
        i, obj(i), delta, size(A,2), nnz(S), subIter);
    
    % testing performance
    Time(i) = toc(t);
    if(isfield(para, 'test'))
        tempS = eye(size(U, 2), size(V, 2));
        RMSE(i) = MatCompRMSE(V, U, tempS, ...
            para.test.row, para.test.col, para.test.data);
    end
    
    if(abs(delta) < tol)
        break;
    end
end

S = eye(size(U, 2), size(V, 2));

out.obj = obj(1:i);
out.rank = nnz(S);
out.RMSE = RMSE(1:i);
out.Time = Time(1:i);

end

